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Abstract : Accuracy of a simulation is strongly depend on the grid quality. Here, quality 
means orthogonality at the boundaries and quasi-orthogonality within the critical regions, 
smoothness, bounded aspect ratios, solution adaptive behaviour, etc. We review various 
functionals for generating high quality structured quadrilateral meshes in two dimensional 
domains. Analysis of Winslow and Modified Liao functionals are presented. Numerical 
>• | examples are also presented to support our theoretical analysis. We will demonstrate the 

use of the Area functional for generating adaptive quadrilateral meshes. 
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1 Introduction 

J> \ Accuracy of numerical solutions of partial differential equations on a grid is very much 

depend on the quality of the underlying grid. There are various parameters for measuring 
grid quality. For example, orthogonality of grid lines and grid density in the regions of 
large solution gradients. A desired grid may be an orthogonal grid with high grid density 
in the areas of sharp solution gradients. Variational methods has been used for improving 
quality of a given grid 0]. In the variational methods, a grid functional is defined. Grid 
functional is an algebraic expression of the position vectors of the internal nodes of a mesh. 
Optimization of the grid functional may result in a grid with desired properties such as 
orthogonal grid lines, equal cell areas, linear or parallelogram cells [see |2| and untangled 
mesh MMJy|- There are many algebraic functionals for grid generation and optimization 
[cf. ^.JSIaHSI^- The first study of grid generation by algebraic functionals were done 
in [lpj. Castillo and Steinberg introduced Length, Orthogonality and Area functionals 
[Tol . Area functional are well known for producing robust quadrilateral meshes. For a 
detailed description of properties of area functionals, the interested readers are referred 



J ■= [gl g2 

9 = J' J 




v*^0) Do) 



gl 



Figure 1: Quantities of interest for 
a quadrilateral cell. 



Figure 2: 2D Structured Mesh. 
Node k is surrounded by four 
quadrilaterals. 



to |9[. Recently the area functional has been used for generating adaptive quadrilateral 
meshes Il2l. 



Let rf) and y(£, rj) be the coordinates of a node in a mesh. Let us further assume 
that x and y are twice differentiate functions of the independent variables £ and rj. An 
integral functional X can be defined as follows 

z(x,y) d = V, x , V, a*. x v , y„) d£ *7 • (1) 

V [0,1] x [0,1] 

We are interested in finding the functions rj) and y(£,rj) for which the integral func- 
tional X attains an extremal value. Such coordinates x and y define a mesh with desirable 
properties. The integral functional X is also referred to as control function for adaptive 
grid generation ij . The conditions for the extremal value of the integral functional X are 
expressed by the Euler-Lagrange equations. The two Euler-Lagrange equations are 



dx d£ drj \dx n 



, (2) 
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dy <9£ V%/ d V \dy v J 

The functions x and y, which satisfy the above Euler-Lagrangian equations, are called 
the extremals of the integral functional X. 

Let us define some quantities of interest. Figure [T] shows a quadrilateral cell, and this 
cell belongs to a mesh. Let the co-variant vector at the node o and in the direction oa 



is gi, and another co- variant vector at the node o but in the direction ob is g2- These 
vectors are given as 

gi = (x a ~ x Q , y a - Vof and g 2 = (x b - x , y b - y Y . (4) 

Other interesting quantities such as the Jacobian and g-tensor matrix can be defined from 
the co-variant vectors. The columns of the Jacobian matrix are the co-variant vectors. 
The g-tensor matrix is the product of the Jacobian matrix with it's transpose. Thus, the 
Jacobian matrix and the g-tensor at the node o and for the cell shown in the Figure are 
given as 

J = [gi g 2 ] and g = J l J. (5) 

The layout of the paper is as follows. In the Section|21 several functionals are presented. 
Continuous and discrete versions of the functionals are presented. Section |3] presents 
several numerical experiments, and finally Section H] concludes the paper. 

2 Discrete Functionals 

Let us first introduce some quantities of interest. These will be used later in formulating 
algebraic functionals. Figure El is a 2 x 2 structured mesh. We use this figure for defining 
these quantities. 

J(ki) refers to the Jacobian (determinant of the Jacobian matrix J (hi)) at the node k 
and for the cell i. Table Q]lists the Jacobian matrix for the four cells surrounding the node 
k- gi(fet) refers to the co- variant base vector at the node k and for the cell i. The base 
vector gx points along horizontal grid lines. Similarly, g2(fc«) refers to the co-variant base 
vector at the node k and for the cell i, and it points along the vertical grid lines. Table 
El lists the co- variant vectors for the Figure 121 It should be noted that column vectors of 
the Jacobian matrix are the co-variant base vectors. For example, the column vectors of 
J(fci) are gi(fci) and g a (fci). That is J(fci) = [gi(fci) ga(fci)]- 

g(ki) refers to the co- variant metric tensor at the node k and for the cell i. It is defined 
as g(ki) = J(k i ) t J(ki). g mn (ki) refers to the (m,n) coefficient of the co-variant metric 
tensor g(ki) for the node k and for the cell i. It can be seen that gn(ki) = gi(fcj)* • gi(fcj) 
and guih) = gi(fcj)* • g2(fcj). Similarly, other coefficients can be defined. 

The coefficient is a measure of the angle between the co-variant base vectors gi 
and g 2 . While, the coefficient gn is a measure of the discrete L 2 length of the co- variant 
vector gi. 

Let us consider a structured quadrilateral mesh (each internal node is surrounded 
by four quadrilaterals) consisting of n internal nodes. The following functionals can be 
defined 



Table 1: Jacobian matrix at the node k for the surrounding cells for the Figure El 
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Table 2: Co- variant vectors at the node k for the surrounding cells for the Figure 121 
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2.1 Area Functional 

The integral form of the standard Area functional is given as 



X A = f i / \J\ 2 d^drj , (6) 

1 J [0,1] x [0,1] 

= / (x^y v - x v y^)dCdr] . (7) 

-/ ro.1ixro.11 



' [0,1] x [0,1 

The Euler-Lagrangian equations for the Area functional are 

^(\J\x v )-^-(\J\x,) = , (8) 

^(\J\y v )--^(\J\yt) = o ■ (9) 

In the simplified form the above equations can be written as 

y 2 x# - x v y v y# - 2.0 j/ { y v x^ + (x e y v + x v y ( ) y^ + y^ 2 x m - x ( y$ y m = , (10) 
x v 2 y# - x v y v x K - 2.0 x^ x v y^ + (x e y v + x r , y$) x (v + x ( 2 y m - y e y s x m = , (11) 

[seelsf. The above Euler-Lagrangian equations are non-elliptic, coupled and quasi-linear 
[cf. |9|. For generating adaptive mesh, the author proposed the following variation in the 



Area functional 



^A(x,y) = 



k=l 



J2 s ^ [ J (**)] s 



(12) 



In the above equation, s(k) is called the adaptive function, and s(h) is the value of 
the adaptive function at the node k and for cell i. 



2.2 Length Functional 

The integral form of the Length functional is given as 

def 1 



X, 



.9u +922} d£dr] , 
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(13) 
(14) 



[(x^) 2 + (x v ) 2 + (y^) 2 + (y v ) 2 ] d^dr] , 

[0,1] x [0,1] 

I; |i| 0, and references therein]. The conditions of extremality of the above length func- 
tional are given by the following Euler-Lagrangian equations 



d 2 x d 2 x 

Q£2 Qyj2 
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Q£2 Q^2 



(15) 
(16) 



The above Laplace's equations can be solved in the computational domain [0, 1] x [0, 1] 
with a specified value of x and y on the boundary. The Euler-Lagrangian equations 
associated with the Length functional are linear and uncoupled. 
The discrete Length functional [l| [l(j is give as follows 



J- L (x,y) = ^ 
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(17) 



2.3 Orthogonality Functional 

The integral form of the Orthogonality functional 0; S S an d references therein] is given 



as follows 
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(19) 
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The Euler-Lagrangian equations corresponding to the minimization of the above integral 
are 
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(21) 
(22) 



These Euler-Lagrangian equations are quasilinear, coupled and non-elliptic in nature 
A simplified form the above Euler-Lagrangian equations is 

x 2 x^ + x v y v y# + (4 x 5 x v + 2 y^ y v ) + y n + x v y^) y^ + x 2 x m + x^ y% y m = 
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(23) 



(24) 



[see |2( • This functional takes only non-negative values, and it would attain a minimum 
value of zero for a completely orthogonal grid. The discrete version of the above Orthog- 
onality functional 0; 0] is given as follows 



^o(x,y) = 
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(25) 



It is found [cf. Q; 0; S 113 3 li | that a linear combination of Area, Length and Orthog- 
onality functionals can produce robust grids in complicated 2D domains. 



2.4 Combination of Length, Area and Orthogonality Function- 
als 

A combined functional is given as 

.F(x,y) =k A J 7 A(x,y) +k h Fi(x,y) +k ^b(x,y) , (26) 

0; 13 II 03 El Q- Here, the parameters k A , k L and k Q satisfy: k A + k L + k G = 1.0 
and k A > 0, kL > 0, ko > 0. A serious drawback of the above combined functional is a 
suitable choice of the parameters. It requires an experience in coming up with a good set 
of parameters n|. It was found k| I3| that the following choice of parameters 



k A = 0.50, k L = 0.0, and k Q = 0.50, (27) 

produces robust grid in many practical domains. The corresponding functional is referred 
as the Knupp's functional. Presented numerical work shows that this functional can 
produce good grids. The Euler-Lagrangian ^13 equations for the minimization of the 



Knupp's functional are 



(x 2 + y 2 ) x^ + 4 x € x v x iv + 2(x 5 y v + x n y e ) y iv + (x 2 + y^ 2 ) x m = , (28) 
(V + ?A, 2 ) y« + 4 y € y„ y € „ + 2(x 5 ^ + x v y{) x in + (x 2 + ^ 2 ) y w = . (29) 



2.5 Winslow Functional 

The Winslow functional is given as follows 



^(x,y) = £ 
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ffll(fcj) + 922(h 
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(30) 



Q;SE3- Here, |J(^)| is the determinant of the Jacobian matrix. One very important 
feature of the above functional is that it has barrier. It means the value of the functional 
approaches infinity when the cells degenerate. That is | J\ — > 0. Thus, this functional 
produces unfolded grids. Numerical experiments also prove this feature of the Winlow 
functional. Since g\\ = gi • gi, and #22 = g2 • g2- It can be shown that the numerator 
(gn + (722) in the Winslow functional (J5U|) is the Frobenius norm of the Jacobian matrix. 
That is 

2 2 

giAh) + g 22 {h) = ^^(J m „(^)) 2 = (|| J(kOII) 2 , 

n=l m=l 

Here, J mn are the components of the Jacobian matrix J. Thus, the Winslow functional 
(|5Ujl can be written as follows 



•F(x,y) = £ 
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IjW 



(31) 



It can be seen easily that the Frobenius norm a 2 x 2 matrix A, and its inverse are 
||^4|| 

related as ||A. _1 || = . Here, \A\ is the determinant of the matrix A. The condition 

number K(A) of a matrix A can be written as K(A) = ||A|| Here, the norm is 

the Frobenius norm. Thus, the Winslow functional can be written as follows 



•F(x,y) = £ 
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(32) 



Thus, the minimization of the functional (|3T)|) is equivalent to the minimization of the 
condition number of the Jacobian matrix. A detailed description of the above analysis 
can also be found in |H H; 13 EiJ- The condition number K(J r (fcj)) can also be expressed 
as 



K(J(fci)) 



\gi(k) x ga(fei)| 



The g(ki) tensor matrix is give as 

g(h) 



gn(h) 912(h) 
921(h) 922(h) 



Let Ai and A2 be the eigenvalues of the matrix g(h). Then 

9n(h) + 922(h) Ai + A 2 
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> 2.0 . 



Here, we have used the relation |«7| 2 = \g\. Thus, the Winslow functional is bounded from 
below. 



2.6 Liao Functional 



The Liao functional for grid generation was proposed in and is give as follows 
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(33) 



2.7 Modified Liao Functional 

The Liao functional can produce folded grids. We will explore it through numerical 
experiments. In the literature, following modification |j| of the Liao functional is given 
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In the above equation, g(ki) is the determinant of the covariant metric tensor g(h)- It 
can be shown that g = J 2 , where J is the Jacobian (determinant of the Jacobian matrix), 
and g is the determinant of the co-variant metric tensor. Thus, this functional, similar 
to the Winslow functional (|3Up. has a barrier. The value of the functional approaches 
infinity when the cells degenerate. That is \J\ — > 0. Thus, this functional produces 
unfolded grids. Numerical experiments also prove this feature of the functional. The 
above functional can remove the folded grids produced by the Liao functional. The 
Modified Liao functional can also be written as follows 



•F(x,y) = ^ 
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(35) 



The Modified Liao functional minimizes the square of the condition number where as the 
Winslow functional minimizes the condition number. 



Figure 3: Example (|3.1|) : Adapted 
Grid by Area Functional. Adapted 
Functional is give as s(x, y) = 5.0 + 
200.0 | sin(2^x) sin(2 7ry)|. 



Figure 4: Example (|3.1|) : Adapted 
Grid by Area Functional. Adapted 
Functional is give as s(x, y) = 5.0 + 
200.0 [sin(2 7rx) wn(2iry)]. 



3 Numerical Examples 

We are interested in finding such a mesh for which the gradient of the functionals vanish. 
The minimization of functionals can be performed by a line search algorithms such as 
Newton's iteration. For the numerical experiments instead of performing the global opti- 
mization we solved the local minimization problems for a single node at a time 0]. In 
all numerical examples, initial grid was generated by linear transfmite interpolation []]]. 

3.1 Adaptive Grid by Area Functional 

It is generally not recommended to uniformly refine the whole mesh in the hope of cap- 
turing the underlying physics. It is desired to adapt a given grid to the requirement of 
the underlying problem. A grid generating algorithm should be able to allocate more grid 
nodes in the part of the domain where large solution gradients occur, and fewer grid nodes 
in the part of the domain where solution is flat. Such grids are called solution- adaptive. 
Behaviour of the underlying solution can be obtained by posteriori indicators • These 
indicators can be computed on a coarse mesh, and they can be used to formulate adaptive 
function s(x,y) in the equation (fT2|) . In the present work, the adaptive function s(x,y) 
is given in the analytic form. 

Figures El and |U report the outcome of our numerical experiments. It should be noted 
that even after adaptation the quadrilateral meshes are convex. One other advantage of 
mesh adaptation by Area functional is that it preserves the mesh topology, and writing a 
solver for a structured mesh is easier compared to unstructured mesh. 



Figure 5: Example (|3.2|) : Folded Grid by Figure 6: Example ()3.2|) : Smooth Grid by 
Transfinite Interpolation. Winslow Functional. 



3.2 Winslow Functional vs Algebraic Method 

Algebraic grid generation methods such as transfinite interpolations 0] are extensively 
used for generating grids. Though, they are one of most simplest method of grid generation 
but algebraic methods can produce folded grids for curved domains as can be seen in 
the Figure 03 One other disadvantage of algebraic grid generation is that boundary 
discontinuity can prorogate inside the domain. It is clear from Figure |H1 that Winslow 
functional smooth the grid, and removes the folded grid lines. 



3.3 Liao, Modified Liao and Area Functionals 

In this example, we perform experiments for comparing Liao, Modified and Area func- 
tional on a simple domain. Outcome of our results are shown in Figures 0IH1 and El It can 
be seen from these figures that Modified Liao functional does indeed removes the inverted 
elements from the mesh but still the quality of the mesh generated by the area functional 
shown in the Figure 01 is certainly better than both Liao and Modified Liao. 



3.4 Length, Area and Knupp's Functionals 

In this example, we compare the Length, the Area and the Knupp functionals. Figures ITTH 
[Til and El are the outcome of our numerical work. The Figure El is a grid by the Length 
functional, the Figure ^2 is a grid by the Area functional, and the Figure ^] is a grid by 
the Knupp's functional. It can be seen that grid by the Area and Knupp's functional are 
better than the grid produced by the Length functional. 



Figure 7: Example ()3.3|) : 
Folded Grid by the Liao 
Functional. 



Figure 8: Example (|3.3|) : 
Grid by the Modified Liao 
Functional. 



Figure 9: Example ()3.3|) : 
Grid by the Area Functional 
with s(x, y) = 1.0. 






Figure 10: Example (|3~4|) : 
Grid by the Length func- 
tional. 



Figure 11: Example (J3~4|) : 
Grid by the Area Functional 
with s(x, y) = 1.0. 



Figure 12: Example (J3~4j) : 
Grid by the Knupp's func- 
tional. 



Figure 13: Example ()3.5j) : Grid by the Figure 14: Example (|3.5j) : Grid by the 
Length Functional. Knupp's Functional. 



3.5 Length and Knupp's Functionals 

In this example, we are comparing Length, and the Knupp's Functional. Outcome of our 
numerical work is reported in Figures EH and Figure EH is a grid generated by the 
Length functional. Figure ITU is grid generated by the Knupp's functional. It can be seen 
that the grid generated by the Knupp's functional is of superior quality. 



4 Conclusions 

We have presented the formulation of various functionals for generating quadrilateral 
meshes, and an analysis of Winslow and Modified Liao functionals that is consistent with 
the numerical experiments. Numerical experiments show that Winslow and Modified 
Liao functionals can remove the folded grids as was expected from theoretical analysis. 
It has been shown that Area functionals can be used for generating robust adaptive 
meshes. Further research is required in formulating adaptive function from a posteriori 
error estimators. 



References 

[1] J.F. Thompson, B.K. Soni and N.P. Weatherill. Handbook of Grid Generation. CRC 
Press, 1998. 

[2] S.K. Khattri. A New Smoothing Algorithm for Quadrilateral and Hexahedral Meshes, 
Lecture Notes in Computer Science, Volume 3992, Apr 2006, Pages 239 - 246, DOI 
10.1007/11758525_32, URL [http://dx.doi.org/10.1007/11758525_32] 



[3] A. M. Winslow. Equipotential zoning of two dimensional meshes. J. Computat. 
Physics, Vol. 1, 1967. 

[4] P.M. Knupp. Jacobian- Weighted Elliptic Grid Generation. SIAM J. Sci. Comput., 
Vol. 17, No. 6, 1475-1490, 1996. 

[5] P.M. Knupp. Algebraic mesh quality metrics. SIAM J. Sci. Comput., Vol. 23, 193- 
218, 2001. 

[6] P.M. Knupp. Hexahedral Mesh Untangling and Algebraic Mesh Quality Metrics. 
Proceedings, 9th International Meshing Roundtable, Sandia National Laboratories. 
Vol. 9, 173-183, 2002. 

[7] J.G. Tinoco-Ruiz and P. Barrer a- Sanchez. Area control in generating smooth and 
convex grids over general plane regions. J. Comput. Appl. Math., Vol. 103, 19-32, 
1999. 

[8] J.G. Tinoco-Ruiz and P. Barrer a- Sanchez. Smooth and convex grid generation over 
general plane regions. Math. Comput. Simulation, Vol. 46, 87-102, 1998. 

[9] J.G. Tinoco, P. Barrera and A. Cortes. Some Properties of Area Functionals in 
Numerical Grid Generation. X Meshing RoundTable, Newport Beach, California, 
USA., 2001. 

[10] J.E. Castillo. On Variational Grid Generation. Ph. D. Thesis, The University of New 
Mexico, Albuquerque, New Mexico, 1987. 

[11] G. Liao and H. Liu. Existence and C 1 ^^ regularity of a minimum of a functional 
related to grid generation problems. Num. Math. PDEs., Vol. 9, 1993. 

[12] S.K. Khattri. An Effective Quadrilateral Mesh Adaptation. Submitted. Available at 
http:/ /www. mi. uib.no/~sanjay/publicatins. html 

[13] P. M. Knupp and S. Steinberg. Fundamentals of Grid Generation. CRC Press. Boca 
Raton, Ann Arbor, London, Tokyo, 1993. 

[14] P. M. Knupp. A robust elliptic grid generator. J. Comp. Phys., 100, 409-418, 1992. 

[15] S. K. Khattri. Hexahedral mesh by area functional. [CA] Simos, Theodore S. (ed.) 
et al., ICNAAM 2005. International conference on numerical analysis and applied 
mathematics 2005. Official conference of the European Society of Computational 
Methods in Sciences and Engineering (ESCMSE), Rhodes, Greek, September 16-20, 
2005. Weinheim: Wiley- VCH. 309-313 (2005). [ISBN 3-527-40652-2/hbk] 

[16] P.M. Knupp, L. Margolin, and M. Shashkov. Reference Jacobian 

Optimization-Based Rezone Strategies for Arbitrary Lagrangian Eule- 

rian Methods. Report No. LA-UR-01-1194. On line avaibale report at 
"http:/ /cnls. lanl.gov/ shashkov/papers/report.ps" . 



S. K. Khattri. Numerical Analysis of an Adaptive Finite Volume Method for 
Single Phase Flow in Highly Heterogenous Medium. Submitted. Available at 
http:/ /www. mi. uib.no/~sanjay/publicatins. html 



o. 

Grid Generation and Adaptation by Functionals 

Sanjay Kumar Khattri 

t^. ■ Department of Mathematics, The University of Bergen, Norway 

sanj ay@mi . uib . no 
http : // www . mi . uib . no/ ~san j ay 



< 



Abstract : Accuracy of a simulation is strongly depend on the grid quality. Here, quality 
means orthogonality at the boundaries and quasi-orthogonality within the critical regions, 
smoothness, bounded aspect ratios, solution adaptive behaviour, etc. We review various 
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1 Introduction 

J> \ Accuracy of numerical solutions of partial differential equations on a grid is very much 

depend on the quality of the underlying grid. There are various parameters for measuring 
grid quality. For example, orthogonality of grid lines and grid density in the regions of 
large solution gradients. A desired grid may be an orthogonal grid with high grid density 
in the areas of sharp solution gradients. Variational methods has been used for improving 
quality of a given grid 0]. In the variational methods, a grid functional is defined. Grid 
functional is an algebraic expression of the position vectors of the internal nodes of a mesh. 
Optimization of the grid functional may result in a grid with desired properties such as 
orthogonal grid lines, equal cell areas, linear or parallelogram cells [see |2| and untangled 
mesh MMJy|- There are many algebraic functionals for grid generation and optimization 
[cf. ^.JSIaHSI^- The first study of grid generation by algebraic functionals were done 
in [lpj. Castillo and Steinberg introduced Length, Orthogonality and Area functionals 
[Tol . Area functional are well known for producing robust quadrilateral meshes. For a 
detailed description of properties of area functionals, the interested readers are referred 
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Figure 1: Quantities of interest for 
a quadrilateral cell. 



Figure 2: 2D Structured Mesh. 
Node k is surrounded by four 
quadrilaterals. 



to |9[. Recently the area functional has been used for generating adaptive quadrilateral 
meshes Il2l. 



Let rf) and y(£, rj) be the coordinates of a node in a mesh. Let us further assume 
that x and y are twice differentiate functions of the independent variables £ and rj. An 
integral functional X can be defined as follows 

z(x,y) d = V, x , V, a*. x v , y„) d£ *7 • (1) 

V [0,1] x [0,1] 

We are interested in finding the functions rj) and y(£,rj) for which the integral func- 
tional X attains an extremal value. Such coordinates x and y define a mesh with desirable 
properties. The integral functional X is also referred to as control function for adaptive 
grid generation ij . The conditions for the extremal value of the integral functional X are 
expressed by the Euler-Lagrange equations. The two Euler-Lagrange equations are 



dx d£ drj \dx n 



, (2) 
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dy <9£ V%/ d V \dy v J 

The functions x and y, which satisfy the above Euler-Lagrangian equations, are called 
the extremals of the integral functional X. 

Let us define some quantities of interest. Figure [T] shows a quadrilateral cell, and this 
cell belongs to a mesh. Let the co-variant vector at the node o and in the direction oa 



is gi, and another co- variant vector at the node o but in the direction ob is g2- These 
vectors are given as 

gi = (x a ~ x Q , y a - Vof and g 2 = (x b - x , y b - y Y . (4) 

Other interesting quantities such as the Jacobian and g-tensor matrix can be defined from 
the co-variant vectors. The columns of the Jacobian matrix are the co-variant vectors. 
The g-tensor matrix is the product of the Jacobian matrix with it's transpose. Thus, the 
Jacobian matrix and the g-tensor at the node o and for the cell shown in the Figure are 
given as 

J = [gi g 2 ] and g = J l J. (5) 

The layout of the paper is as follows. In the Section|21 several functionals are presented. 
Continuous and discrete versions of the functionals are presented. Section |3] presents 
several numerical experiments, and finally Section H] concludes the paper. 

2 Discrete Functionals 

Let us first introduce some quantities of interest. These will be used later in formulating 
algebraic functionals. Figure El is a 2 x 2 structured mesh. We use this figure for defining 
these quantities. 

J(ki) refers to the Jacobian (determinant of the Jacobian matrix J (hi)) at the node k 
and for the cell i. Table Q]lists the Jacobian matrix for the four cells surrounding the node 
k- gi(fet) refers to the co- variant base vector at the node k and for the cell i. The base 
vector gx points along horizontal grid lines. Similarly, g2(fc«) refers to the co-variant base 
vector at the node k and for the cell i, and it points along the vertical grid lines. Table 
El lists the co- variant vectors for the Figure 121 It should be noted that column vectors of 
the Jacobian matrix are the co-variant base vectors. For example, the column vectors of 
J(fci) are gi(fci) and g a (fci). That is J(fci) = [gi(fci) ga(fci)]- 

g(ki) refers to the co- variant metric tensor at the node k and for the cell i. It is defined 
as g(ki) = J(k i ) t J(ki). g mn (ki) refers to the (m,n) coefficient of the co-variant metric 
tensor g(ki) for the node k and for the cell i. It can be seen that gn(ki) = gi(fcj)* • gi(fcj) 
and guih) = gi(fcj)* • g2(fcj). Similarly, other coefficients can be defined. 

The coefficient is a measure of the angle between the co-variant base vectors gi 
and g 2 . While, the coefficient gn is a measure of the discrete L 2 length of the co- variant 
vector gi. 

Let us consider a structured quadrilateral mesh (each internal node is surrounded 
by four quadrilaterals) consisting of n internal nodes. The following functionals can be 
defined 



Table 1: Jacobian matrix at the node k for the surrounding cells for the Figure El 



J(h) 



[X4 


- Xk) 


(Xl 


- x k ) 


(Va 


- Vk) 


(2/1 


-Vk) 


(x 2 


- x k ) 


(x 3 


- x k ) 


(2/2 


- Vk) 


(Vz 


-Vk) 



J(k 2 ) 



(x 2 


- x k ) 


(xi 


- x k ) 




-Vk) 


(yi 


- Vk) 


(x 4 


- x k ) 


(x 3 


- x k ) 


(2/4 


-Vk) 


(2/3 


- Vk) 



Table 2: Co- variant vectors at the node k for the surrounding cells for the Figure 121 



giifa) 





- x k 


V\ 


- Vk 


x 2 


- x k 


2/2 


- Vk 



12(h) 

12(h) 



Xl - x k 
V\ ~ Vk 

X3 - x k 
2/3 - Vk 



5i(fe) 
li(h) 



x 2 


- x k 


2/2 


- Vk 


X4 


- x k 


2/4 


- Vk 



12(h) 
52(^4) 



Xl - x k 
2/1 - Vk 

X3 - x k 
2/3 - Vk 



2.1 Area Functional 

The integral form of the standard Area functional is given as 



X A = f i / \J\ 2 d^drj , (6) 

1 J [0,1] x [0,1] 

= / (x^y v - x v y^)dCdr] . (7) 

-/ ro.1ixro.11 



' [0,1] x [0,1 

The Euler-Lagrangian equations for the Area functional are 

^(\J\x v )-^-(\J\x,) = , (8) 

^(\J\y v )--^(\J\yt) = o ■ (9) 

In the simplified form the above equations can be written as 

y 2 x# - x v y v y# - 2.0 j/ { y v x^ + (x e y v + x v y ( ) y^ + y^ 2 x m - x ( y$ y m = , (10) 
x v 2 y# - x v y v x K - 2.0 x^ x v y^ + (x e y v + x r , y$) x (v + x ( 2 y m - y e y s x m = , (11) 

[seelsf. The above Euler-Lagrangian equations are non-elliptic, coupled and quasi-linear 
[cf. |9|. For generating adaptive mesh, the author proposed the following variation in the 



Area functional 



^A(x,y) = 



k=l 



J2 s ^ [ J (**)] s 



(12) 



In the above equation, s(k) is called the adaptive function, and s(h) is the value of 
the adaptive function at the node k and for cell i. 



2.2 Length Functional 

The integral form of the Length functional is given as 

def 1 



X, 



.9u +922} d£dr] , 



[0,1] x [0,1] 



(13) 
(14) 



[(x^) 2 + (x v ) 2 + (y^) 2 + (y v ) 2 ] d^dr] , 

[0,1] x [0,1] 

I; |i| 0, and references therein]. The conditions of extremality of the above length func- 
tional are given by the following Euler-Lagrangian equations 



d 2 x d 2 x 

Q£2 Qyj2 

d 2 y d 2 y 

— - H = . 

Q£2 Q^2 



(15) 
(16) 



The above Laplace's equations can be solved in the computational domain [0, 1] x [0, 1] 
with a specified value of x and y on the boundary. The Euler-Lagrangian equations 
associated with the Length functional are linear and uncoupled. 
The discrete Length functional [l| [l(j is give as follows 



J- L (x,y) = ^ 



k=l 



^2 {gnih) + 922(h)) 



i=l 



(17) 



2.3 Orthogonality Functional 

The integral form of the Orthogonality functional 0; S S an d references therein] is given 



as follows 



J, 



def 



O 



(g 12 ) 2 d£drj 



[0,1] x [0,1] 



,1 ■ 62J 



: d£ drj, 



[0,1] x [0,1] 



[0,1] x [0,1] 



(x ( x v + y^y v ) 2 d^dr] 



(18) 
(19) 
(20) 



The Euler-Lagrangian equations corresponding to the minimization of the above integral 
are 



d 



dx\ d 



ot, \ or] J or) \ ot, J 







9l2 



dy 



d£ V 9rj J dr] 



d ( dy 

+ — I 9\2 



9i2^r = , 







(21) 
(22) 



These Euler-Lagrangian equations are quasilinear, coupled and non-elliptic in nature 
A simplified form the above Euler-Lagrangian equations is 

x 2 x^ + x v y v y# + (4 x 5 x v + 2 y^ y v ) + y n + x v y^) y^ + x 2 x m + x^ y% y m = 



y v 2 y& + Vv x v x tt + ( 4 Vi Vv + 2 x ? ^) %»? + ^ + ^ x c) %7 + 2/? 2 ftm + % x ? x 



'/'/ 



(23) 



(24) 



[see |2( • This functional takes only non-negative values, and it would attain a minimum 
value of zero for a completely orthogonal grid. The discrete version of the above Orthog- 
onality functional 0; 0] is given as follows 



^o(x,y) = 



k=l 



i=i 



(25) 



It is found [cf. Q; 0; S 113 3 li | that a linear combination of Area, Length and Orthog- 
onality functionals can produce robust grids in complicated 2D domains. 



2.4 Combination of Length, Area and Orthogonality Function- 
als 

A combined functional is given as 

.F(x,y) =k A J 7 A(x,y) +k h Fi(x,y) +k ^b(x,y) , (26) 

0; 13 II 03 El Q- Here, the parameters k A , k L and k Q satisfy: k A + k L + k G = 1.0 
and k A > 0, kL > 0, ko > 0. A serious drawback of the above combined functional is a 
suitable choice of the parameters. It requires an experience in coming up with a good set 
of parameters n|. It was found k| I3| that the following choice of parameters 



k A = 0.50, k L = 0.0, and k Q = 0.50, (27) 

produces robust grid in many practical domains. The corresponding functional is referred 
as the Knupp's functional. Presented numerical work shows that this functional can 
produce good grids. The Euler-Lagrangian ^13 equations for the minimization of the 



Knupp's functional are 



(x 2 + y 2 ) x^ + 4 x € x v x iv + 2(x 5 y v + x n y e ) y iv + (x 2 + y^ 2 ) x m = , (28) 
(V + ?A, 2 ) y« + 4 y € y„ y € „ + 2(x 5 ^ + x v y{) x in + (x 2 + ^ 2 ) y w = . (29) 



2.5 Winslow Functional 

The Winslow functional is given as follows 



^(x,y) = £ 



fc=i 



4 

E 

1=1 



ffll(fcj) + 922(h 

\J(h)\ 



(30) 



Q;SE3- Here, |J(^)| is the determinant of the Jacobian matrix. One very important 
feature of the above functional is that it has barrier. It means the value of the functional 
approaches infinity when the cells degenerate. That is | J\ — > 0. Thus, this functional 
produces unfolded grids. Numerical experiments also prove this feature of the Winlow 
functional. Since g\\ = gi • gi, and #22 = g2 • g2- It can be shown that the numerator 
(gn + (722) in the Winslow functional (J5U|) is the Frobenius norm of the Jacobian matrix. 
That is 

2 2 

giAh) + g 22 {h) = ^^(J m „(^)) 2 = (|| J(kOII) 2 , 

n=l m=l 

Here, J mn are the components of the Jacobian matrix J. Thus, the Winslow functional 
(|5Ujl can be written as follows 



•F(x,y) = £ 



k=l 



E 

i=l 



IjW 



(31) 



It can be seen easily that the Frobenius norm a 2 x 2 matrix A, and its inverse are 
||^4|| 

related as ||A. _1 || = . Here, \A\ is the determinant of the matrix A. The condition 

number K(A) of a matrix A can be written as K(A) = ||A|| Here, the norm is 

the Frobenius norm. Thus, the Winslow functional can be written as follows 



•F(x,y) = £ 



k=l 



i=l 



(32) 



Thus, the minimization of the functional (|3T)|) is equivalent to the minimization of the 
condition number of the Jacobian matrix. A detailed description of the above analysis 
can also be found in |H H; 13 EiJ- The condition number K(J r (fcj)) can also be expressed 
as 



K(J(fci)) 



\gi(k) x ga(fei)| 



The g(ki) tensor matrix is give as 

g(h) 



gn(h) 912(h) 
921(h) 922(h) 



Let Ai and A2 be the eigenvalues of the matrix g(h). Then 

9n(h) + 922(h) Ai + A 2 



\J(h)\ 



V Ai Aj 



> 2.0 . 



Here, we have used the relation |«7| 2 = \g\. Thus, the Winslow functional is bounded from 
below. 



2.6 Liao Functional 



The Liao functional for grid generation was proposed in and is give as follows 



k=l 



J2^n + 922 2 + 2 g 12 2 



i=i 



(33) 



2.7 Modified Liao Functional 

The Liao functional can produce folded grids. We will explore it through numerical 
experiments. In the literature, following modification |j| of the Liao functional is given 



•F(x,y) = £ 



fe=i 



i=l 



9n(h) + 922(h 

\/g(h) 



(34) 



In the above equation, g(ki) is the determinant of the covariant metric tensor g(h)- It 
can be shown that g = J 2 , where J is the Jacobian (determinant of the Jacobian matrix), 
and g is the determinant of the co-variant metric tensor. Thus, this functional, similar 
to the Winslow functional (|3Up. has a barrier. The value of the functional approaches 
infinity when the cells degenerate. That is \J\ — > 0. Thus, this functional produces 
unfolded grids. Numerical experiments also prove this feature of the functional. The 
above functional can remove the folded grids produced by the Liao functional. The 
Modified Liao functional can also be written as follows 



•F(x,y) = ^ 



k=l 



i=i 



(35) 



The Modified Liao functional minimizes the square of the condition number where as the 
Winslow functional minimizes the condition number. 



Figure 3: Example (|3.1|) : Adapted 
Grid by Area Functional. Adapted 
Functional is give as s(x, y) = 5.0 + 
200.0 | sin(2^x) sin(2 7ry)|. 



Figure 4: Example (|3.1|) : Adapted 
Grid by Area Functional. Adapted 
Functional is give as s(x, y) = 5.0 + 
200.0 [sin(2 7rx) wn(2iry)]. 



3 Numerical Examples 

We are interested in finding such a mesh for which the gradient of the functionals vanish. 
The minimization of functionals can be performed by a line search algorithms such as 
Newton's iteration. For the numerical experiments instead of performing the global opti- 
mization we solved the local minimization problems for a single node at a time 0]. In 
all numerical examples, initial grid was generated by linear transfmite interpolation []]]. 

3.1 Adaptive Grid by Area Functional 

It is generally not recommended to uniformly refine the whole mesh in the hope of cap- 
turing the underlying physics. It is desired to adapt a given grid to the requirement of 
the underlying problem. A grid generating algorithm should be able to allocate more grid 
nodes in the part of the domain where large solution gradients occur, and fewer grid nodes 
in the part of the domain where solution is flat. Such grids are called solution- adaptive. 
Behaviour of the underlying solution can be obtained by posteriori indicators • These 
indicators can be computed on a coarse mesh, and they can be used to formulate adaptive 
function s(x,y) in the equation (fT2|) . In the present work, the adaptive function s(x,y) 
is given in the analytic form. 

Figures El and |U report the outcome of our numerical experiments. It should be noted 
that even after adaptation the quadrilateral meshes are convex. One other advantage of 
mesh adaptation by Area functional is that it preserves the mesh topology, and writing a 
solver for a structured mesh is easier compared to unstructured mesh. 



Figure 5: Example (|3.2|) : Folded Grid by Figure 6: Example ()3.2|) : Smooth Grid by 
Transfinite Interpolation. Winslow Functional. 



3.2 Winslow Functional vs Algebraic Method 

Algebraic grid generation methods such as transfinite interpolations 0] are extensively 
used for generating grids. Though, they are one of most simplest method of grid generation 
but algebraic methods can produce folded grids for curved domains as can be seen in 
the Figure 03 One other disadvantage of algebraic grid generation is that boundary 
discontinuity can prorogate inside the domain. It is clear from Figure |H1 that Winslow 
functional smooth the grid, and removes the folded grid lines. 



3.3 Liao, Modified Liao and Area Functionals 

In this example, we perform experiments for comparing Liao, Modified and Area func- 
tional on a simple domain. Outcome of our results are shown in Figures 0IH1 and El It can 
be seen from these figures that Modified Liao functional does indeed removes the inverted 
elements from the mesh but still the quality of the mesh generated by the area functional 
shown in the Figure 01 is certainly better than both Liao and Modified Liao. 



3.4 Length, Area and Knupp's Functionals 

In this example, we compare the Length, the Area and the Knupp functionals. Figures ITTH 
[Til and El are the outcome of our numerical work. The Figure El is a grid by the Length 
functional, the Figure ^2 is a grid by the Area functional, and the Figure ^] is a grid by 
the Knupp's functional. It can be seen that grid by the Area and Knupp's functional are 
better than the grid produced by the Length functional. 



Figure 7: Example ()3.3|) : 
Folded Grid by the Liao 
Functional. 



Figure 8: Example (|3.3|) : 
Grid by the Modified Liao 
Functional. 



Figure 9: Example ()3.3|) : 
Grid by the Area Functional 
with s(x, y) = 1.0. 






Figure 10: Example (|3~4|) : 
Grid by the Length func- 
tional. 



Figure 11: Example (J3~4|) : 
Grid by the Area Functional 
with s(x, y) = 1.0. 



Figure 12: Example (J3~4j) : 
Grid by the Knupp's func- 
tional. 



Figure 13: Example ()3.5j) : Grid by the Figure 14: Example (|3.5j) : Grid by the 
Length Functional. Knupp's Functional. 



3.5 Length and Knupp's Functionals 

In this example, we are comparing Length, and the Knupp's Functional. Outcome of our 
numerical work is reported in Figures EH and Figure EH is a grid generated by the 
Length functional. Figure ITU is grid generated by the Knupp's functional. It can be seen 
that the grid generated by the Knupp's functional is of superior quality. 



4 Conclusions 

We have presented the formulation of various functionals for generating quadrilateral 
meshes, and an analysis of Winslow and Modified Liao functionals that is consistent with 
the numerical experiments. Numerical experiments show that Winslow and Modified 
Liao functionals can remove the folded grids as was expected from theoretical analysis. 
It has been shown that Area functionals can be used for generating robust adaptive 
meshes. Further research is required in formulating adaptive function from a posteriori 
error estimators. 
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